Numerical studies of Casimir interactions 
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We study numerically the Casimir interaction between dielectrics in both two and three dimensions. We 
demonstrate how sparse matrix factorizations enable one to study torsional interactions in three dimensions. 
In two dimensions we study the full cross-over between non-retarded and retarded interactions as a function 
of separation. We use constrained factorizations in order to measure the interaction of a particle with a rough 
dielectric surface and compare with a scaling argument. 
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Dispersion forces have as their origin the fluctuations of po- 
larization in materials coupled to the long-ranged electrody- 
namic interaction described by Maxwell's equations. The first 
convincing demonstration of their importance was the calcu- 
lation by Kees5)m of the interaction between fluctuating clas- 
sical dipoles ill]. The introduction of quantum fluctuations by 
London accounted for the long-ranged, part of the 

van der Waals interaction in most materials. Later, Casimir 
and Polder ijsfl showed that retardation modifies the interac- 
tions in an important manner- leading to a decay in the in- 
teraction which is asymptotically l/?-^ at zero temperature. 
Further advances were made by the Russian school 0] who 
showed how to formulate the interactions in terms of the di- 
electric response of materials. Overviews with many refer- 
ences to theoretical and experimental developments are to be 
found in isl S 01 ■ Retarded Casimir interactions are the dom- 
inant interaction between neutral surfaces at the submicron 
scale. 

Whilst the analytic basis of the theory is largely estab- 
lished its application is difficult in experimentally interest- 
ing geometries. One is constrained to work withperturba- 
tive expansion about exactly solvable geometries [81], or use 
ad hoc schemes such as the proximity force approximation. 
Only a few geometries have been attacked with exact analytic 
techniques 191]. Recently several attempts have been made 
to study numerically the interactions by using methods from 
modern com puta tional science- including fast multigrid lat- 
tice solvers fllOll in order to calculate Green functions and 
forces, or the use of discretized determinants to determine free 
energies ll 111 . 

In this Letter we will present a series of techniques which 
enable one to evaluate the interaction between dielectric bod- 
ies in full vectorial electrodynamics. Firstly, we calculate the 
torsional potential between two three-dimensional bodies in 
the retarded regime, using a full discretization of Maxwell's 
equations, we note that the Casimir torque has recently re- 
ceived the attention of experimentalists Il2lll3ll . For more de- 
tailed studies we present results for two-dimensional systems. 
This allows us to study the cross-over between the near- and 
far-field regimes and also to measure the interaction between 
a particle and a rough surface. With these two-dimensional 
systems we implement general strategies which substantially 
increase the efficiency of simulations, at the same time de- 
creasing the sensitivity of the results to numerical round-off 
errors. 



In three dimensions we discretize Maxwell's equations to a 
cubic Yee-lattice ifTill . lattice constant a = 1, associating the 
electric degrees of freedom to the links, magnetic degrees of 
freedom are localized on the faces of the lattice. We remind 
the reader that the finite difference approximation to the V x 
operator, here designated Curl, maps the electric field on four 
links surrounding the face of the cube to the magnetic field. 
Curl is needed in the Maxwell equation 
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= -cCurlE 



The adjoint operator maps fields from the faces to the links. 
We will denote it Curl*. It intervenes in the second time de- 
pendent Maxwell equation. 
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The importance of clearly distinguishing the two operators 
will become apparent when we discuss the two-dimensional 
case below. We use Heaviside-Lorentz units in which 
Maxwell's equations are directly parameterized by the speed 
of light in vacuum, c. 



From these two equations Lifshitz theory H15I1 shows that 
the free energy of interaction between dielectric bodies is 
found from from the imaginary time wave equation for the 
vector potential in the temporal gauge where E = — A/c and 



Cm-rCurU A = Pa A = 



e(r, 



Alternatively one introduces a magnetic formulation and 
works with a potential such that H = G/c and considers the 
wave equation 



Curl 
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e(r,w) 



Curl* \ G = VaG = Q 



In our work we always consider the differences in free en- 
ergy between pairs of configurations; we thus avoid a full ac- 
count of the self-energy variations of dielectric media ll 1 1. 
The free energy difference between two configurations 1 , 2 is 
found from 



j ^ {Indet pi (w) -Indet (i) 



FIG. 1: A Pair of structured dielectric rings. Eacli quadrant has 
different dielectric properties. 
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for either choice of wave operator, Va or Vq', while self- 
energy contributions are different in the two formulations we 
have verified with our codes that both give the same result for 
the long-ranged part of the interactions that we are interested 
in. 

We perform the frequency integration in eq. ([T]) by chang- 
ing integration variables to z, where uj = az/{\ — z) with 
< z < 1. The parameter a is chosen so that the ma- 
jor features in the integrand occur for values of z near 1/2. 
We then use iVg -point Legendre-Gauss quadrature to replace 
the integral by a weighted sum over discrete frequencies. 
We evaluate determinants by finding the Cholesky factoriza- 
tion Ld of 'D{lo) such that ijj is lower triangular lfl6ll and 
LdLJ) ~ 'Di'jj)- The determinant of T) is then given by 

hidetPH = 2^1n(LD,,,,) 

i 

When we examine the detailed structure of Maxwell's 
equations discretized lo V = sites in three dimensions 
we discover that the Curl operator is a matrix of dimension 
?>V X 3y and has 12V non-zero elements. The operator 
(Curl* Curl) has 39^ non-zero elements. The major tech- 
nical difficulty comes from the fact that the matrices we work 
with have dimensions which are very large, ^ 10^ x 10^. 
All numerical work was performed with an Intel Xeon-5140 
workstation. 

We now calculate the Casimir torque between two parallel 
rings centered on a common axis, figure [T] Each ring is di- 
vided into quadrants with alternating dielectric properties. We 
take permittivities which are independent of frequency, corre- 
sponding to the full retarded regime 11511 with ei(Lj) = 5, 
e2{^) — 10; the space around the rings is a vacuum with 

= 1- We measure the energy of interaction as the top ring 
is rotated with respect to the lower. The zero of the inter- 
action corresponds to aligned rings. As the rings are rotated 
the interface between the dielectric materials, as interpolated 
to the lattice, undergoes some re-arrangement changing the 
self energy of the rings. We thus perform two runs. The first 
run of a single rotating ring determines this variation in the 
self-energy. The second run with the both rings allows one to 
measure the interaction energy by subtraction. 

We worked with a system of dimensions y = 55 x 55 x 55, 
figure 121 The graph of the interaction energy as a function of 
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FIG. 2: Interaction energy as a function of angle as a ring of figure[T] 
rotates. In tlie linear parts of tlie curve tlie torque is almost indepen- 
dent of the angle. Ring diameter 36a, separation and thickness 2a. 
Rounding is determined by the ratio of separation to diameter of the 
rings. 10 days of computation with Ng = 8. Cholesky factor 9GB. 



angle is noticeably triangular in shape between 7r/8 to Btt/S. 
This is understood by the fact that the interaction energy is 
dominated by the interactions directly across the gap. The 
fluctuations in the curve about the expected linear behavior, 
together with its slight asymmetry give an idea of the noise 
coming from irregularities of the interpolation of the disks to 
the lattice. This irregularity is particularly clear on comparing 
the points for 7r/4 and 37r/4. 

We now turn to two-dimensional electrodynamics where we 
can study systems with larger linear dimensions. Such large 
system sizes are needed in order to follow the cross-overs be- 
tween different regimes in the interaction of particles or if one 
wishes to simulate structured or disordered materials in order 
understand the efficiency of analytic approximations. 

In three dimensions the two formulation in terms of Va 
and Vg are largely equivalent. In two-dimensional electrody- 
namics this is no longer the case. Consider an electrodynamic 
system in which there are two components of the electric field 
in the x — y plane; the magnetic field then has just a single 
component in the z direction. The Curl operator becomes a 
rectangular matrix of dimensions yx2y where now V ~ L^. 
The standard formulation in terms of the vector potential gives 
to an operator Va of dimensions 2V x 2V with lAV non-zero 
elements; the alternative formulation in terms of Vq leads to 
determinants of dimensions V x V involving just 5V non- 
zero elements; the size of the matrix that we must work with 
is smaller in the 2?g formulation. We used Dq in the follow- 
ing numerical work, having checked that we obtain equivalent 
results. 

We started by measuring the cross-over between the short- 
ranged non-retarded interaction to the long-ranged Casimir 
force. We studied a pair of dielectric particles described by 



3 



the single pole approximation to the dielectric constant 

e(w) = 1 + ^ , „. „ 

where x is the zero frequency electric susceptibility. The in- 
teraction is retarded for separations D ^ c/ojq, non-retarded 
for D <^ c/uoo. 

We measured the interaction between two dielectric par- 
ticles in a square, periodic cell of dimensions L x L using 
SuiteSparse llTIl to perform both the ordering and the factor- 
ization of the matrices. We placed a first particle at the ori- 
gin, and considered two possible positions of a second parti- 
cle to calculate a free energy difference using eq. ([T]i. The first 
results were disappointing- rather small systems (L = 50) 
were sensitive to numerical round-off errors. The origin of 
this problem was quite clear. In a large system there is an ex- 
tensive self-energy ^ L^. Pair interactions calculated as the 
difference between two large numbers are unreliable. 

We avoided this problem by separating the free energy con- 
tributions from the neighborhood of the three interesting sites 
and the rest of the system. We did this by introducing a 
block-wise factorization of V that enabled us to both solve 
the round-off problem while re-using much of the numerical 
effort need to generate the Cholesky factors thus improving 
the efficiency of the code. 

We now write the symmetric matrix from the wave equa- 
tion in block form, V = (^y'^ ' determinant 

is dct(2?) = dct(X) det(S') where the Schur complement 
S = Z - Y'^X^^Y 111]. We group sites so that the great 
majority is within the block X and sites that we are interested 
in are in the block Z. It is the term in dct{X) the gives the 
large extensive free energy which caused our numerical prob- 
lems. It is independent of the properties of our test particles. 
All the interesting information on energy differences is in the 
Schur complement, S. 

We start by finding the Cholesky factorization of X, L^- 
The Schur complement is calculated by solving the triangular 
equations LxU = Yhy forward substitution, then calculating 
S = Z — U'^U. Our separation of energies into an extensive 
constant and a small set of interacting sites allows us to study 
the interaction of systems of sizes up to L = 2000 before 
round-off becomes a problem. 

In order to generate data we generalized the method to a 
three level scheme- firstly collect the set of sites (here ~ 100) 
of all the separations required to generate a curve into the 
block Z, and form the Schur complement forming a small ef- 
fective theory for all these remaining sites. Within the smaller 
matrix that has been generated we again re-order to succes- 
sively put each interesting sets of variables in the bottom-right 
corner of the effective theory and find the Schur complement 
of these remaining variables. We can then calculate interac- 
tions between the particles while minimizing round-off errors. 

We remind the reader that in two dimensions the electro- 
static potential is logarithmic between two charges, and that 
dipole-dipole fluctuations lead to van der Waals interactions 
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FIG. 3: Scaled interaction free energy, —Ur'^ for a pair of dielectric 
particles (e(0) = 8) in a box of dimensions 2000 x 2000 as a function 
of separation. Curves from top to bottom correspond to luq/c = 
10, 0.3, 0.1, 0.03, 0.01 0.003. For large loq/c, Ur^ is constant, □. 
For smaller coo /ewe see both retarded and non-retarded interactions. 
Solid line corresponds to U^dw ~ l/r*. 10GB for Cholesky factor. 
Six hours of calculation. Ng = 25. 

decaying as Uydw = As in three dimensions retardation 
leads to an accelerated decay so that the Casimir interaction 
varies as Uc ^ 1/r^. In our simulations we used values of 
luq/c varying from 0.003 to 10, figure[3] We determined the 
energy of interaction of particles U, as a function of separa- 
tion r while moving the second particle in the simulation cell 
out to {L/5, L/5); the zero of energy is calculated for two 
particles separated by {L/2, L/2). We scale out the retarded 
behavior, plotting —U{r)r^. We see that for the largest ujq/c 
the interactions are retarded for all separations, □. For the 
smaller values of ujq/c the interaction varies as 1 /?■'*. In the 
scaled curve this gives the linear rise clearly visible in the fig- 
ure, . For 0.1 < ujq/c < 0.01 we see both the near- and 
far-field behaviors clearly displayed within a single sample- 
permitting the detailed study of cross-over phenomena with 
frequency dependent dielectric behavior No assumptions of 
symmetry are made in the calculation; the method can be used 
with bodies of arbitrary geometry. 

We now turn to a problem where analytic results are much 
more difficult to find: The interaction of a dielectric particle 
with a rough surface, figured We generated rough surfaces as 
realizations of solid-on-solid random walks on a lattice. Ap- 
proximately half of the simulation box contains dielectric ma- 
terial with e = 8,^0 oo; the rest of the box has e = 1. 
We measure the interaction with a test particle as a function 
of the distance from the rough surface using the above method 
of block Schur complements to perform a single large factor- 
ization per frequency for each realization of the disorder We 
generated 1000 rough surfaces and measured the average in- 
teraction with the surface (U), as a function of separation, as 
well as the variance in the potentials. 

We understand the results, figure |5] with a scaling argu- 
ment. When the particle is a distance r from the surface the 
interaction is dominated by a front of length r along the sur- 
face. Since the surface is a random walk its average posi- 
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FIG. 4: Realization of rough interface and set of measurement posi- 
tions, X , for tlie interaction energy whicli will be separated into the 
block Z. Anisotropic horizontal and vertical scales. 



ble with m. 

We have demonstrated the power of direct methods from 
Unear algebra when applied to the study of dispersion forces. 
In three dimensions we have measured interactions in experi- 
mentally realizable geometries- though system sizes are still 
too small to accurately measure cross-overs between differ- 
ent scaling regimes. In two dimensions we have shown how 
to measure the cross-over between London dispersion and 
Casimir interactions, and have determined correction to scal- 
ing exponents for the interactions of a disordered systems. 

Work financed in part by Volkswagenstiftung. 




FIG. 5: (I) o, —{U)r^, averaged interaction between dielectric par- 
ticle and rough dielectric surface. (2) o, —UsT'^, interaction between 
particle and flat surface. (3) A, (7„r'^, variance of interaction for 
rough surfaces. (4) U, SUr^, difference in mean interaction en- 
ergy between a flat and a rough surface. Solid lines: r~"^'^ and r~*. 
L — 1000. Two weeks of simulation time. Cholesky factor 2.5GB. 

Ng = 20. 

tion is displaced by Sr ^ ±r^/^ compared to the flat surface. 
The interaction between a smooth surface and a particle varies 
as Us ^ l/r^ in the Casimir regime. The interaction of the 
particle should thus he U ^ l/{r + 6r)^. If we expand to 
first order we find that the variance of the interaction should 
scale as, (A) (t„ ^ ,,-3.5 yyjjjjg ^jje second order expansion 
gives a shift in the mean potential, ([/), which varies as, (□), 
6U ^ l/r**. The numerical data are compatible with this scal- 
ing. The argument is easily generalized to affine surfaces with 
other, less trivial roughness exponents giving results compati- 
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